function [status_quo] = func_allocations_imputed( IJipref,IJipref_choice,IJoutput, IJclassSize,IJVAm1, IJVAm2,IJfracDisadv, IJitidx,IJjtidx,IJimputedteacher,IJappyear,yearselect,makebalanced,firstbest)

% This computes teacher propose stable allocations and first best 
% This code clears the market in two rounds:
% (1) among the non-imputed teachers
% (2) among any unfilled positions, it clears it with imputed teachers

% makebalance=1/0 specifices whether we make imputed teachers balanced with
% positions, or whether we make imputed teachers long 
% (note that non-imputed teachers are thus always short)

% Output is: (1) output for disadvantaged; (2) output for advantaged
% (3) output for disadvantaged among imputed teachers; (4) output for
% advantaged among imputed teachers

% The reason to look at output among imputed teachers is that making those
% teachers balanced or long potentially affects things

%if firstbest==1
%status_quo=[va_disadv_fb va_adv_fb va_disadv_imp_fb va_adv_imp_fb]; 

%end 

%if firstbest==0
%status_quo=[va_disadv va_adv va_disadv_imp va_adv_imp ]; 

%end 



   IJnonimputed=logical(1-IJimputedteacher);

    %% go down to IJappyear==yearselect
    IJipref=IJipref(IJappyear==yearselect);
    IJipref_choice=IJipref_choice(IJappyear==yearselect);
    IJitidx2015=IJitidx(IJappyear==yearselect);
    IJjtidx2015=IJjtidx(IJappyear==yearselect);
    IJclassSize=IJclassSize(IJappyear==yearselect);
    IJoutput=IJoutput(IJappyear==yearselect);
    IJimputedteacher=IJimputedteacher(IJappyear==yearselect);
    IJnonimputed=IJnonimputed(IJappyear==yearselect);
    
    IJVAm1=IJVAm1(IJappyear==yearselect);
    IJVAm2=IJVAm2(IJappyear==yearselect);
    IJfracDisadv=IJfracDisadv(IJappyear==yearselect);

    
    % redo indexes so that they start at 1 
    
    IJitidx=IJitidx2015-min(IJitidx2015)+1;
    IJjtidx=IJjtidx2015-min(IJjtidx2015)+1;
    
    
      for i=1:length(IJitidx)
           IJnonimputed_sq(IJitidx(i,1),IJjtidx(i,1))=IJnonimputed(i,1);    % these are I X J 
      end

    
    if makebalanced==1  % here, makebalanced is dropping teachers (i); want to keep the non-imputed teachers, and only drop imputed teachers 
        jmax=max(IJjtidx);       
        ilist=unique(IJitidx);
        nonimputelist=ilist(IJnonimputed_sq(:,1));
        imputelist=ilist(logical(1-IJnonimputed_sq(:,1)));
        jmax_sub=jmax-length(nonimputelist); 
        ikeep_impute=randsample(imputelist,jmax_sub);  % randomly sample imputed teachers
        ikeep=[nonimputelist; ikeep_impute]; 
        iselect=ismember(IJitidx,ikeep);  % generate a logical vector telling us which to keep
        
        % now apply the selector
        
        IJipref=IJipref( iselect);
        IJipref_choice=IJipref_choice( iselect);       
        IJclassSize=IJclassSize( iselect);
        IJoutput=IJoutput( iselect);
        IJVAm1=IJVAm1( iselect);
        IJVAm2=IJVAm2( iselect);
        IJfracDisadv= IJfracDisadv( iselect);
        IJitidx=IJitidx( iselect);
        IJjtidx= IJjtidx( iselect);    
        IJimputedteacher=IJimputedteacher(iselect);
        
        
        % sort i in sequential order (to renumber), then reorder everything
        % else
        [b1,Jsorter]=sort(IJitidx); % sort by i
        inbetween=repmat((1:jmax),jmax,1);
        IJitidx=inbetween(:);
        
        IJipref=IJipref(Jsorter);
        IJipref_choice=IJipref_choice(Jsorter);        
        IJclassSize=IJclassSize( Jsorter);
        IJoutput=IJoutput( Jsorter);
        IJVAm1=IJVAm1( Jsorter);
        IJVAm2=IJVAm2(Jsorter);
        IJfracDisadv= IJfracDisadv(Jsorter);
        IJjtidx=IJjtidx(Jsorter);   
        IJimputedteacher=IJimputedteacher(Jsorter);  

    end
    
  clearvars IJnonimputed_sq     
    
    
%% make the teacher preferences rectangular, and then rank





   for i=1:length(IJitidx)
         IJjpref_sq(IJjtidx(i,1),IJitidx(i,1))=IJoutput(i,1);    % these are J X I (useful for DA)
         IJfracDisadv_ij(IJjtidx(i,1),IJitidx(i,1))=IJfracDisadv(i,1);  % these are J x 1 (useful for first best)   
         IJclassSize_ij(IJjtidx(i,1),IJitidx(i,1))=IJclassSize(i,1);  % these are I x J (useful for first best)
         IJijpref_sq(IJitidx(i,1),IJjtidx(i,1))=IJoutput(i,1);    % these are I X J (useful for first best)
         IJclassSize_sq(IJitidx(i,1),IJjtidx(i,1))=IJclassSize(i,1);  % these are I x J (useful for first best)
         IJVAm1_sq(IJitidx(i,1),IJjtidx(i,1))=IJVAm1(i,1);  % these are I X J (useful for first best)
         IJVAm2_sq(IJitidx(i,1),IJjtidx(i,1))=IJVAm2(i,1);  % these are I X J (useful for first best)
         IJfracDisadv_sq(IJitidx(i,1),IJjtidx(i,1))=IJfracDisadv(i,1);  % these are I X J (useful for first best)  
         IJimputedteacher_sq(IJitidx(i,1),IJjtidx(i,1))=IJimputedteacher(i,1);  % these are I X J (useful for first best)  
         
         IJipref_sq(IJitidx(i,1),IJjtidx(i,1))=IJipref(i,1); % I x J 
         IJipref_choice_sq(IJitidx(i,1),IJjtidx(i,1))=IJipref_choice(i,1); %   I X J          
            
   end
   
   IJtype1num=(1-IJfracDisadv_sq).*IJclassSize_sq;
   IJtype2num=IJfracDisadv_sq.*IJclassSize_sq;
   
   % split teacher preferences into stuff for imputed and non-imputed
   % teachers
   
   IJimputed=logical(IJimputedteacher_sq(:,1));
   IJnonimputed=logical(1-IJimputed);

   
if firstbest==0    
   IJipref_choice_sq_nonimp=IJipref_choice_sq(IJnonimputed,:);
   
   % split position preferences into stuff for nonimputed
   
   IJjpref_sq_nonimp=IJjpref_sq(:,IJnonimputed'); % note: this is J X I 
      
    % do this for nonimputed teachers
    [B1, IJirank_nonimp]=sort(IJipref_choice_sq_nonimp,2,'descend');  % this is sorted in descending order; so the first column in I is the most preferred of the teacher in the relevant row 

   % do this over nonimputed teachers 
    
    [B1, IJjrank_nonimp]=sort(IJjpref_sq_nonimp,2,'descend');  % this is sorted in descending order; so the first column in I is the most preferred of the school in the relevant row 
    
    % make the school output vector rectangular, and then rank
    % note: this breaks ties in an arbitrary way 
   
   % run DA
          
   [jass_iprop_nonimp]= galeshapley_gender_unbalanced(IJirank_nonimp,IJjrank_nonimp);  % teach propose
   
   jass_iprop_nonimp_wide=[ jass_iprop_nonimp, (1:length(jass_iprop_nonimp))']; % teacher index, then school index (i,j)   

    % identify unmatched positions as those where the teacher index is 0,
    % and create the matrix of matches to add later
    
    unmatched=(jass_iprop_nonimp_wide(:,1)==0);   
    matched=logical(1-unmatched);
    unmatched_school=jass_iprop_nonimp_wide(unmatched,2);
    jass_iprop_nonimp_wide_final=jass_iprop_nonimp_wide(matched,:);
    
    
    % construct output stuff for the nonimputed teachers 
 
    
    IJVAm2_sq_nonimp=IJVAm2_sq(logical(IJnonimputed),:);
    IJVAm1_sq_nonimp=IJVAm1_sq(logical(IJnonimputed),:);
    IJtype1num_nonimp=IJtype1num(logical(IJnonimputed),:);  
    IJtype2num_nonimp=IJtype2num(logical(IJnonimputed),:);   

    
   for i=1:length(jass_iprop_nonimp_wide_final)
          i_type2_tot_nonimp(i,1)=IJVAm2_sq_nonimp(jass_iprop_nonimp_wide_final(i,1), jass_iprop_nonimp_wide_final(i,2)).* IJtype2num_nonimp(jass_iprop_nonimp_wide_final(i,1), jass_iprop_nonimp_wide_final(i,2));
          i_type2_stu_nonimp(i,1)= IJtype2num_nonimp(jass_iprop_nonimp_wide_final(i,1), jass_iprop_nonimp_wide_final(i,2));
          
          i_type1_tot_nonimp(i,1)=IJVAm1_sq_nonimp(jass_iprop_nonimp_wide_final(i,1), jass_iprop_nonimp_wide_final(i,2)).* IJtype1num_nonimp(jass_iprop_nonimp_wide_final(i,1), jass_iprop_nonimp_wide_final(i,2));
          i_type1_stu_nonimp(i,1)= IJtype1num_nonimp(jass_iprop_nonimp_wide_final(i,1), jass_iprop_nonimp_wide_final(i,2));
   end
   
    va_disadv_nonimp=mean(IJVAm2_sq_nonimp(jass_iprop_nonimp_wide_final(:,1),1));
    va_adv_nonimp=mean(IJVAm1_sq_nonimp(jass_iprop_nonimp_wide_final(:,1), 1));
    %va_disadv_nonimp=sum(i_type2_tot_nonimp)/sum(i_type2_stu_nonimp);
    %va_adv_nonimp=sum(i_type1_tot_nonimp)/sum(i_type1_stu_nonimp);
   
   % create versions of these matrices for the unmatched folks: these are
   % imputed teachers (IJimputed) and unmatched positions (unmatched)
   
   IJVAm2_sq_open=IJVAm2_sq(logical(IJimputed),unmatched');
   IJVAm1_sq_open=IJVAm1_sq(logical(IJimputed),unmatched');
   IJtype2num_open=IJtype2num(logical(IJimputed),unmatched');
   IJtype1num_open=IJtype1num(logical(IJimputed),unmatched');
    
    % now create preferences among the imputed teachers for the open
    % positions, and similarly preferences among open positions for the
    % imputed teachers
    
    IJipref_choice_sq_imp=IJipref_choice_sq(logical(IJimputed),unmatched');
    IJjpref_sq_imp=IJjpref_sq(unmatched,logical(IJimputed)'); % note: this is J X I 
   
   [B1, IJirank_imp]=sort(IJipref_choice_sq_imp,2,'descend');  % this is sorted in descending order; so the first column in I is the most preferred of the teacher in the relevant row     
   [B1, IJjrank_imp]=sort(IJjpref_sq_imp,2,'descend');  % this is sorted in descending order; so the first column in I is the most preferred of the school in the relevant row 
 
   % run DA
   
   [jass_iprop_imp]= galeshapley_gender_unbalanced(IJirank_imp,IJjrank_imp);  % teach propose
   jass_iprop_imp_wide=[ jass_iprop_imp, (1:length(jass_iprop_imp))']; % teacher index, then school index (i,j)   
   
   % now compute output among the imputed teachers
   
   
   
   
   for i=1:length(jass_iprop_imp_wide)
          i_type2_tot_imp(i,1)=IJVAm2_sq_open(jass_iprop_imp_wide(i,1), jass_iprop_imp_wide(i,2)).* IJtype2num_open(jass_iprop_imp_wide(i,1), jass_iprop_imp_wide(i,2));
          i_type2_stu_imp(i,1)= IJtype2num_open(jass_iprop_imp_wide(i,1), jass_iprop_imp_wide(i,2));
          
          i_type1_tot_imp(i,1)=IJVAm1_sq_open(jass_iprop_imp_wide(i,1), jass_iprop_imp_wide(i,2)).* IJtype1num_open(jass_iprop_imp_wide(i,1), jass_iprop_imp_wide(i,2));
          i_type1_stu_imp(i,1)= IJtype1num_open(jass_iprop_imp_wide(i,1), jass_iprop_imp_wide(i,2));
   end
   
   
   % given that the marginal matched people are the imputed, it is of
   % interest to compute (1) overall output and (2) output among imputed
   % (and see how that varies with market balance)
   
   % among imputed
   va_disadv_imp=mean(IJVAm2_sq_open(jass_iprop_imp_wide(:,1),1));
   va_adv_imp=mean(IJVAm1_sq_open(jass_iprop_imp_wide(:,1), 1));
 
   %va_disadv_imp=sum(i_type2_tot_imp)/sum(i_type2_stu_imp);
    %va_adv_imp=sum(i_type1_tot_imp)/sum(i_type1_stu_imp);
    
    % total
    
    i_type2_tot=[i_type2_tot_nonimp;  i_type2_tot_imp];
    i_type2_stu=[i_type2_stu_nonimp;  i_type2_stu_imp];
    
    i_type1_tot=[i_type1_tot_nonimp;  i_type1_tot_imp];
    i_type1_stu=[i_type1_stu_nonimp;  i_type1_stu_imp];

    va_disadv=sum(i_type2_tot)/sum(i_type2_stu);
    va_adv=sum(i_type1_tot)/sum(i_type1_stu);

end           
  % compute first-best maximizing disadvantaged output
  
if firstbest==1
    
    
    % first do first-best among nonimputed teachers
    
    IJijpref_sq_nonimp=IJijpref_sq(logical(IJnonimputed),:);
    
    out_sq=-IJijpref_sq_nonimp;
    [iind_best,jind_best] = linear_sum_assignment(out_sq);
        
    jass_iprop_wide=[ iind_best,jind_best]; % teacher index, then school index (i,j)   
  
   % figure out unassigned positions
   
   jmax=max(IJjtidx); 
   schoollist=(1:jmax)';
   unmatched=setdiff(schoollist,jind_best);
   
   IJijpref_sq_open=IJijpref_sq(logical(IJimputed),unmatched');
   
   
    IJVAm2_sq_nonimp=IJVAm2_sq(logical(IJnonimputed),:);
    IJVAm1_sq_nonimp=IJVAm1_sq(logical(IJnonimputed),:);
    IJtype1num_nonimp=IJtype1num(logical(IJnonimputed),:);  
    IJtype2num_nonimp=IJtype2num(logical(IJnonimputed),:); 
    
    
      %  output among nonimputed
   for i=1:length(jass_iprop_wide)
          i_type2_tot_nonimp(i,1)=IJVAm2_sq_nonimp(jass_iprop_wide(i,1), jass_iprop_wide(i,2)).* IJtype2num_nonimp(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
          i_type2_stu_nonimp(i,1)= IJtype2num_nonimp(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
          
          i_type1_tot_nonimp(i,1)=IJVAm1_sq_nonimp(jass_iprop_wide(i,1), jass_iprop_wide(i,2)).* IJtype1num_nonimp(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
          i_type1_stu_nonimp(i,1)= IJtype1num_nonimp(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
   end
   
    % among nonimputed
    va_disadv_nonimp_fb=mean(IJVAm2_sq_nonimp(jass_iprop_wide(:,1), 1));
    va_adv_nonimp_fb=mean(IJVAm1_sq_nonimp(jass_iprop_wide(:,1), 1));

   % va_disadv_nonimp_fb=sum(i_type2_tot_nonimp)/sum(i_type2_stu_nonimp);
   % va_adv_nonimp_fb=sum(i_type1_tot_nonimp)/sum(i_type1_stu_nonimp);
   
   % make versions of these input variables among the unmatched
   
   IJVAm2_sq_open=IJVAm2_sq(logical(IJimputed),unmatched');
   IJVAm1_sq_open=IJVAm1_sq(logical(IJimputed),unmatched');
   IJtype2num_open=IJtype2num(logical(IJimputed),unmatched');
   IJtype1num_open=IJtype1num(logical(IJimputed),unmatched');
   
   % now solve for first best among the imputed and remaining positions
   
    out_sq=-IJijpref_sq_open;
    [iind_best,jind_best] = linear_sum_assignment(out_sq);
        
    jass_iprop_wide=[ iind_best,jind_best]; % teacher index, then school index (i,j)    
    
    % now solve for output among the imputed teachers 
 
      for i=1:length(jass_iprop_wide)
          i_type2_tot_imp(i,1)=IJVAm2_sq_open(jass_iprop_wide(i,1), jass_iprop_wide(i,2)).* IJtype2num_open(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
          i_type2_stu_imp(i,1)= IJtype2num_open(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
          
          i_type1_tot_imp(i,1)=IJVAm1_sq_open(jass_iprop_wide(i,1), jass_iprop_wide(i,2)).* IJtype1num_open(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
          i_type1_stu_imp(i,1)= IJtype1num_open(jass_iprop_wide(i,1), jass_iprop_wide(i,2));
      end
   
   
         % among imputed
    va_disadv_imp_fb=mean(IJVAm2_sq_open(jass_iprop_wide(:,1), 1));     
    va_adv_imp_fb=mean(IJVAm1_sq_open(jass_iprop_wide(:,1), 1));     

    %va_disadv_imp_fb=sum(i_type2_tot_imp)/sum(i_type2_stu_imp);
    %va_adv_imp_fb=sum(i_type1_tot_imp)/sum(i_type1_stu_imp);
    
    % total
    
    i_type2_tot=[i_type2_tot_nonimp;  i_type2_tot_imp];
    i_type2_stu=[i_type2_stu_nonimp;  i_type2_stu_imp];
    
    i_type1_tot=[i_type1_tot_nonimp;  i_type1_tot_imp];
    i_type1_stu=[i_type1_stu_nonimp;  i_type1_stu_imp];

    va_disadv_fb=sum(i_type2_tot)/sum(i_type2_stu);
    va_adv_fb=sum(i_type1_tot)/sum(i_type1_stu);
  

end     




%% collect output

if firstbest==1
status_quo=[va_disadv_fb va_adv_fb va_disadv_imp_fb va_adv_imp_fb va_disadv_nonimp_fb va_adv_nonimp_fb]; 

end 

if firstbest==0
status_quo=[va_disadv va_adv va_disadv_imp va_adv_imp va_disadv_nonimp va_adv_nonimp ]; 

end 


end

